Linear scaling computation of the Fock matrix. VIII. 
Periodic boundaries for exact exchange at the F-point. 
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A translationally invariant formulation of the Hartree-Fock (HF) F-point approximation is pre- 
sented. This formulation is achieved through introduction of the Minimum Image Convention (MIC) 
at the level of primitive two-electron integrals, and implemented in a periodic version of the ONX 
algorithm [J. Chem. Phys, 106 9708 (1997)] for linear scaling computation of the exchange ma- 
trix. Convergence of the HF-MIC F-point model to the HF k-space limit is demonstrated for fully 
periodic magnesium oxide, ice and diamond. Computation of the diamond lattice constant using 
the HF-MIC model together with the hybrid PBEO density functional [Theochem, 493 145 (1999)] 
yields ao — 3.569A with the 6-21G* basis set and a 3 x 3 x 3 supercell. Linear scaling computation 
of the HF-MIC exchange matrix is demonstrated for diamond and ice in the condensed phase. 
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I. INTRODUCTION 

In a preceding articlei, methods were introduced 
for constructing the periodic F-point Coulomb and 
Exchange-Correlation matrix that achieve a cost scaling 
only linearly with system size, N. In this paper, an 0{N) 
algorithm is presented for computation of the periodic 
Hartree-Fock (HF) exchange matrix at the F-point. The 
Hartree-Fock approximation is often a fast, first approx- 
imation and also a starting point for correlated "wave- 
function" methods. Also, hybrid Hartree-Fock/Density 
Functional Theory (HF/DFT) model chemistries are an 
important next step in accuracy beyond the Generalized 
Gradient ApproximationSiSiii^. Together with linear scal- 
ing methods for computing the density matrijsSiij these 
advances provide a framework for the application of both 
HF and HF/DFT models to large condensed phase sys- 
tems, surfaces and wires. 

To date, condensed phase HF and HF/DFT calcula- 
tions have been carried out almost exclusively with the 
ab initio solid state program CRYSTAL - , which employs 
conventional algorithms for evaluation of the HF ex- 
change matrix. The CRYSTAL program employs meth- 
ods for k-space integration, which generates wavefunc- 
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tions with the correct translational symmetries. While 
the F-point approximation forgoes k-space integration, 
it should recover the correct symmetries in the limit of a 
large super cell. With a more tractable formulation and 
asymptotic correctness in the limit of large systems, the 
F-point approximation would seem an ideal basis for lin- 
ear scaling exchange algorithms. However in preliminary 
studies, we found that the naive HF F-point approxi- 
mation converges to position dependent values, and that 
even in the large super cell limit, values different from 
the k-space integration limit are approached. 



In this paper we develop a translationally invariant 
definition of F-point Hartree-Fock exchange, which cor- 
rectly approaches the k-space integration value in the 
limit of a large super-cell. This is accomplished through 
introducing a Minimum Image Convention (MIC) into 
the exchange kernel at the level of primitive two-electron 
integrals, as described in Section IIIII Incorporation of 
the MIC integrals into a periodic version of the ONX 
algorithm^ for linear scaling computation of the HF ex- 
change matrix is then described in Section In Section 
IVll we demonstrate convergence of the HF-MIC model 
chemistry to the k-space integration limit for a number 
of systems, and in Section lVHI we demonstrate linear scal- 
ing for both diamond and ice. In Scction FVIIII we discuss 
these results, and then in Section llXI we summarize our 
results. 
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II. PERIODIC EXACT EXCHANGE 



III. EXCHANGE AT THE T-POINT 



In the conventional formulation of periodic boundary 
conditions, the Bloch functions 



J2 e''' ^ 

R 



0a(r-R), 



(1) 



are often constructed from non-orthogonal functions local 
to the unit cell. Here, the local function (j>a is a Gaussian- 
Type Atomic Orbital (GTAO) centered on atom A, while 
the sum on R runs over the Bravais lattice vectors defined 
by integer translates of the lattice primitives a, b and c. 
Inversely, the vector k runs over the reciprocal lattice 
vectors. 

To date, rapid computation of the Hartree-Fock ex- 
change interaction demands the analytic evaluation of 
two-electron integrals, which is possible when the local 
basis functions are of Cartesian Gaussian type. Typi- 
cally, these functions have the form 



Mr) = {x-A^,y-iy-Ayr^iz-A,r^e 



-Ca(r-A)2 



(2) 



where the triad {/a,ma,na} sets angular symmetry and 
the exponent Ca is chosen to describe a particular length 
scale. Gaussian basis functions are often contracted to 
approximate atomic eigenfunctions. 

With periodic boundary conditions, the exact Hartree- 
Fock exchange matrix igSiiii 



K[k]=^K[L]e'''•^ 



where 



(3) 



(4) 



HN,cd 



is defined in real space, and the two-electron integrals, 
written in chemist's notation, are 



drdr' 



, 0a(r)0c(r + H)(/)b(r' + L)Mr' + H + N) 



(5) 



The formally infinite sums over lattice vectors in 
Eq. 10} involve many contributions that are in practice 
infinitesimal. In part this is due to the decay between 
local basis function products (paipc] the product of two 
Guassians centered at A and C decays also as a Gaus- 
sian with separation | A — C| . Truncation based purely on 
the overlap of Gaussian basis functions is reliable and well 
controlled, leading to 0{N) product terms. Additionally, 
the density matrix is known to decay exponentially for 
non-metallic systems. A cutoff defining the range of al- 
lowed exchange interactions may be used to exploit this 
fall off a priori by confining summation over p>jSii2iiiii2ii^ 
to only those terms that satisfy the imposed geometric 
constraints. Methods equivalent to a cutoff have also 
been used to achieve A'^-scaling of the exchange matrix in 
gas phase calculationsi^ii^. 



The F-point approximation limits k-space sampling to 
just the central cell at k = 0. In the naive F-point limit, 
we take P[N] = ^n,o P, reducing Eq. Q to the less cum- 
bersome relation: 



p, 



cd 



aHI 



'^b^d. 



(6) 



HL.crf 



where we have relabeled the lattice sums for clarity. How- 
ever, there are two significant problems with this naive 
approach. First, the exact exchange potential, which is 
periodic in Eq. 01 has been truncated asymmetrically to 
yield Eq. Neighboring cells correctly representing the 
exchange interaction have been removed, leading to a 
non-periodic exchange potential that varies with position 
of the coordinates in relation to the unit cell. Effectively, 
a position dependent environment has been created for 
each charge distribution. Secondly, this truncated ex- 
pression violates symmetry of the generalized exchange 
energy: 



where 



^;,[p^p'']^^;,[p^p^ 



E^[P'',P''] = Tr [P" • K (P*" 



(7) 



(8) 



and the superscripts label different density matrices. 

Let us consider the first problem, which is the most 
serious defect of Eq. |B1 This translational variance 
has a well understood analogue in classical molecular 
dynamicsiSiiiii^, where the arbitrary truncation of short 
range potentials also leads to numerical artifacts in ener- 
gies and forces. In both cases, translational invariance is 
restored by introducing the Minimum Image Gonvention 
(MIG), which ensures that interactions are always cal- 
culated between nearest images. In computation of the 
Hartree-Fock exchange matrix, the MIC F-point approx- 
imation is just 



Kab = --^ E ^"'^ ' 



MIC ■ 



(9) 



HLxd 



where the MIC condition is applied in computation of the 
two-electron integrals at the contraction phase, ensuring 
that primitive charge distributions interact consistently 
over a minimum distance. In particular, if the primitive 
basis function product 0q(/)^ is centered at P and the 
primitive product 4>b'Pd is at Q, then the minimum image 
convention is applied to the interaction vector PQ = 
P — Q using 



pq = PQ 



(10a) 



PQi = PQi- Mim{pqi-SlG}i{5,pqi)) (10b) 
PQmic = M-pq (10c) 

where 6 w 10""'^^ is needed to avoid wrapping errors^^, 
and M is the 3x3 shape matrix of the unit cell, composed 
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of the primitive lattice vectors, 

(ax Cx \ 
ay by Cy (ll) 
ciz bz J 

This approach is completely general, and can be used at 
the primitive level with any modern approach to com- 
puting two-electron integrals. While the MIC yields 
translational invariance of the exchange matrix, it does 
not recover Exchange Kernel Permutational Symmetry 
(EKPS). With greater expense, we could recover the 
EKPS by explicitly symmetrizing the primitive Gaussian 
products within the kernel. However, in the limit of large 
systems, where the range of the density matrix becomes 
smaller than the system size, the EKPS is recovered as 
demonstrated in Section lVll fnote that, without the MIC, 
EKPS is not recovered in this limit). 

IV. OPTIMAL DAMPING AND SYMMETRY 
OF THE EXCHANGE KERNEL 

For difficult, unstable SCF problems, the Optimal 
Damping Algorithm (ODA) of Cancesi^ is an efficient 
method that guarantees convergence of the HE model. 
However, permutational symmetry of the exchange ker- 
nel is an implicit, simplifying assumption in formulation 
of the conventional ODA algorithm. For small periodic 
systems, violation of the EKPS creates problems for the 
ODA, leading to a non-quadratic behavior (the HF model 
should yield an exactly quadratic, convex minimization 
problem). While the EKPS is restored with increasing 
system size, loose numerical thresholds can also lead to 
loss of the EKPS in the limit of a large system. In both 
cases, loss of EKPS can lead to incorrect determination 
of the ODA mixing parameter 

where the superscripts indicate consecutive steps in the 
SCF cycle given by the endpoints, with A G [0,1]. It 
is of course a simple matter to reformulate the ODA, 
using definitions for the endpoint derivatives that do not 
assume EKPS: 

+Tr [piF°] + Tr [P^F^] - 2Tr [P°F°] 

^ ^El~El + Tr[{P^-P^)T] 

-Tr [P^F"] - Tr [P^F^] + 2Tr [P^F^] , 

(13) 

where F is the Fockian, T is the kinetic energy matrix 
and Enc is the nuclear-electrostatic energy. This modified 
ODA leads to the correct quadratic parameterization and 
guarantees convergence of the HF-MIC model. 



V. IMPLEMENTATION 

A general treatment of F-point periodic boundary con- 
ditions has been implemented in the MondoSCF— suite 
of programs for linear scaling quantum chemistry. A 
detailed account of these developments for pure Den- 
sity Functional Theory has been given in a compan- 
ion paper^, including the periodic development of the 
Quantum Chemical Tree Code (QCTC) for A^-scaling 
Coulomb summation. In addition to QCTC, the lin- 
ear scaling, Quartic Trace- ReSetting (TRS4)'' density 
matrix solver has been used throughout, together with 
inverse congruence transformations provided by sparse 
atom-blocked approximate AINV— . 

The Order N eXchange (ONX) algorithm^S for com- 
puting the gas phase exchange matrix has been modified 
by placing dual loops running over the lattice vectors H 
and L around the original ONX loop structures. Two 
ordered bra and ket distribution buffers are assembled 
for each lattice vector pair, which are then used to drive 
the basic ONX algorithm. The Minimum Image Conven- 
tion has been introduced into the primitive contraction 
stage, in the Vertical Recurrence Relations component 
of a symmetry driven Head-Gordon Pople22. scheme for 
computing two-electron integrals. 

All developments were implemented in MoNDoSCF 
v1.0q;9^*', a suite of linear scaling Quantum Chemistry 
codes. The code was compiled using the Portland Group 
F90 compiler PGF90 v4.22^ with the -01 -tp athlon op- 
tions and with the GNU C compiler GCC v3.2.2 using the 
-01 flag. All calculations were carried out on a 1.6GHz 
AMD Athlon running RedHat LiNUXvQ.OSi. 



TABLE I: Progression of Hartree-Fock F-point super-cell cal- 
culations for MgO using the periodic RHF-MIC and RHF 
8-511G/8-51G level of theory. Comparison is made to a fi- 
nal value approaching the k-space integration limit for the 
primitive cell. 



Program 


A^MgO 


Energy (au) 


Energy/ iVMgO 


MondoSCF 




-274.530033 


-137.265017 






-1098.43823 


-137.304778 




16'' 


-2197.05636 


-137.316023 






-4394.48467 


-137.327646 




54' 


-7415.89617 


-137.331411 






-8789.24941 


-137.332022 




128* 


-17578.5022 


-137.332048 




216'= 


-29663.7291 


-137.332079 


CRYSTAL98" 


2d 


-274.664153 


-137.332077 



"F-point 
''Triclinic 
'^Cubic 



■^8 X 8 X 8 k-space integration grid 



TABLE II: Progression of F-point super-cell calculations of 
proton ordered ice at the RHF-MIC/8-51G/5-11G* level of 
theory. Comparison is made to a final value approaching the 
k-space integration limit for the primitive cell. 



Program 


Nh20 


Energy (an) Energy/iVnaO 


DEKPS 


MONDOSCF" 


2 


-152.03025 


-76.01512 


10" 


M 




16 


-1216.3003 


-76.01877 


10" 


4 




54 


-4105.0163 


-76.01882 


10" 


-5 




128 


-9730.4090 


-76.01882 


10" 


-6 




250 


-19004.705 


-76.01882 


10" 


-6 


CRYSTALgS" 


2 


-152.03765 


-76.01882 







"F-point 

'6x6x6 k-space integration grid 



TABLE III: Lattice constants in A for diamond computed 
using different system sizes, theory levels, and basis sets at 
a LOOSE accuracy. For comparison, the experimental value 
for diamond, extrapolated to T = OK, is ao = 3.567A, while 
ao = 3.583A is obtained with both the RPBE/6-31G* GGA 
and the RTPSS/6-31G* meta-GGA level of theory2&. 



Program 




Basis 


ag 




MONDOSCF" 


64 


ST0-3G 


3.587 


3.601 




216 


STO-3G 


3.582 


3.594 


CRYSTAL* 


2 


ST0-3G 


3.581 




MONDOSCF" 


64 


6-21G* 


3.575 


3.571 




216 


6-21G* 


3.571 


3.569 


CRYSTAL'' 


2 


6-21G* 


3.574 





"F-point 
Taken from Ref. 



VI. VALIDATION 

Comparison is made to the Gaussian orbital, periodic 
ab initio program CRYSTAL98^', primarily using basis 
sets optimized for the condensed phase, obtained from 
Ref. 

Table shows the progression of total energies com- 
puted with MONDOSCF for MgO with the F-point RHF- 
MIC model, using the 8-511G basis set for magnesium 
and the 8-51G basis set for oxygen. Comparison is made 
to a final value approaching the k-space integration limit 
for the primitive ceU obtained using CRYSTAL98. The 
CRYSTAL basis sets were obtained from Ref. [2^, and 
the primitive cubic and triclinic cell coordinates used for 
this system are given in Ref. p|. The values controlling 
accuracy of the CRYSTAL98 program were obtained 
from Ref. [sJI, while the MoNDoSCF calculations were 
carried out using the TIGHT level of accuracy, defining 
numerical thresholds delivering numbers precise to the 
digits quoted. 

In Table ^1 total energies computed with the RHF- 
MIC F-point super cell approach are listed for proton 
ordered iceSI, using the 8-51G basis for oxygen and the 
5-llG* basis for hydrogen. A comparison is made to a 



FIG. 1: The relaxed, uniaxial lattice potential of proton or- 
dered ice^'^. Comparison is made between a RHF-MIC/8- 
51G/5-11G* 250 molecule F-point super-cell calculation and 
a CRYSTAL98 calculation carried out with a two molecule 
primitive cell using a 6 x 6 x 6 k-space integration grid. 



^ -76.0188 



Crystal98(6x6x6) 
MondoSCF(250) 




7.1 7.2 

Lattice Constant (A) 



final value obtained using Cyrstal98 with a 6 x 6 x 6 
k-space integration grid as explained in Ref. |23|- The 
MONDOSCF values were obtained using the GOOD accu- 
racy level, defining numerical thresholds that deliver to- 
tal energies precise to the digits quoted. The primitive 
cubic cell coordinates used in these ice calculations are 
listed in Ref. j30] . Also shown in this table is a measure of 
the Deviation from the Exchange Kernel Pcrmutational 
Symmetry (DEKPS), 



DEKPS 



(14) 



where the superscripts refer simply to consecutive steps 
in the SCF cycle. Because this measure is normalized, it 
yields a roughly consistent measure throughout the SCF. 
As shown in Tabled the DEKPS decreases with system 
size to a constant value consistent with the numerical 
threshold used to discard matrix elements. 

Table Hill compares the diamond lattice constant com- 
puted by MondoSCF using the RHF-MIC F-point super- 
cell approximation to results obtained in reference [26|. 
The result achieved using the 216 atom 3x3x3 super- 
cell, the 6-21G* basis set and the hybrid PBEO functional 
yields ao = 3. 569 A and should be compared to the ex- 
perimental, T = OK result of ao = 3.567A, as well as 
the value ao — 3.583A computed with the TPSS meta- 
GGA25. 

Figure^shows the unrelaxed, uniaxial lattice potential 
of proton ordered icc^^ in the a direction (see Ref. [s^l 
for the coordinate system). Comparison is made between 
a 250 molecule RHF-MIC F-point super-cell calculation 
performed with MondoSCF and CRYSTAL98 calcu- 
lations carried out with a two molecule primitive cell 
using a 6 X 6 X 6 k-space integration grid, the 8-51G 
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basis for oxygen and the 5-llG* basis for hydrogen. The 
MONDOSCF GOOD option was used, dehvering a relative 
accuracy of 6-7 digits. The potential minimum for the 
CRYSTAL98 curve is 7.144A and for the MoNDoSCF 
calculation 7. 145 A. 

FIG. 2: CPU time for the exchange, Coulomb and density 
matrix build (scaled by 10) of diamond using LOOSE thresh- 
olding at the RHF-MIC/STO-3G level of theory 



the RHF-MIC model, GOOD thresholds, the 8-51G basis 
set for oxygen and the 5-llG* basis set for hydrogen. 
This is the level of theory used to produce Table m and 
Fig. An early onset of linear scaling is observed at 
about 50 water molecules, which is similar to the onset 
achieved for water clusters in the gas phased. 

VIII. DISCUSSION 



H 
D 

Oh 
U 



□ Exchange Matrix Build 
o Coulomb Matrix Build 
^ Density Matrix Solve (x 10) 




100 200 300 400 500 600 

Number of Carbon Atoms 



800 



FIG. 3: CPU time for the exchange, Coulomb and density 
matrix build (scaled by 10) of proton ordered ice using GOOD 
thresholding at the RHF-MIC/8-51G/5-11G* level of theory 
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Differences between MondoSCF and CRYSTAL98 
include the use of fitting functions and k-space integra- 
tion by CRYSTAL98 and the use of approximate 0{N) 
scaling algorithms and the F-point approximation by 
MondoSCF. Nevertheless, in the large super cell limit, 
we have achieved a very close agreement between the two 
methods. 

Importantly, the periodic MondoSCF algorithms 
have also demonstrated linear scaling for three- 
dimensional systems. While achieving linear scaling in 
three-dimensions for gas phase systems is challenging, it 
is even more difficult with full periodicity. This is espe- 
cially true for dense systems like diamond; linear scaling 
was achieved here only for large unit cells, a minimal ba- 
sis and a LOOSE accuracy level. It is possible to achieve 
an early onset of linear scaling with significantly larger 
Gaussian basis sets by tightening the valence region, but 
doing so without care can lead to poor performance in the 
computation of proprieties such as the lattice constant. 
Localization of the basis set has received recent, criti- 
cal attention^SiSi^, leading to numerical basis functions 
that yield both linear scaling and high quality proper- 
ties. So far though, this issue remains unexplored in the 
context of a Gaussian representation. 

For less dense systems and larger basis functions, the 
periodic MondoSCF algorithms are able to achieve an 
early onset of linear scaling with more accurate thresh- 
olding parameters. With parallelism^^iSSiS, this opens 
the prospect of performing hybrid HF/DFT simulations 
of fluids. 

These differences in behavior highlight the use of con- 
tinuous thresholds that set the cost to accuracy ratio on 
the fly. Thus, as the gap closes or the basis sets becomes 
overcomplete, periodic ONX and TRS4 will correctly re- 
vert to 0{N'^) and 0{N^) algorithms respectively. This 
should be compared with use of a radial cutoff to de- 
termine a priori the graph of the density and exchange 
matrix, which docs not correctly permit fill-in. 



VII. SCALING 

In Fig. 13 timings are shown for building the exact ex- 
change, the Coulomb, and the density matrix of diamond 
at the RHF-MIC/ST0-3G level of theory with a LOOSE 
set of thresholds. This is the same accuracy level used to 
compute the values listed in Table lllll 

Figure |3| shows timings for building the exchange. 
Coulomb and density matrix of proton ordered iceSL at 



IX. CONCLUSIONS 

For the first time, a translationally invariant definition 
of the periodic Hartree-Fock F-point approximation has 
been presented. Based on inserting the Minimum Image 
Condition (MIC) into the contraction phase of periodi- 
cally sumirred two-electron integral algorithms, this HF- 
MIC approach has been used to extend the ONX algo- 
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rithm for linear scaling computation of the exchange ma- 
trix to periodic boundary conditions. Convergence of the 
HF-MIC F-point super-cell approximation to the k-space 
integration limit has been demonstrated for MgO and ice 
to better than 8 digits. Linear scaling was demonstrated 
for diamond and ice, including MIC-exchange, Coulomb 
and density matrix construction. 
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APPENDIX A: THRESHOLDS 

There are four thresholds that control the numeri- 
cal precision of the linear scaling algorithms in Mon- 
doSCF. They are the matrix threshold Tmtrix described 
in Ref. ^3, the two-electron threshold (thresh in 
Ref. M), the distribution threshold Tdist described in 
Ref. 0, and the HiCu threshold Thicu used to set the 
exchange-correlation grid {tt in Ref. i,38j). The thresh- 
old also controls accuracy of the Coulomb matrix as 
discussed in Ref. [l|. These thresholds, listed in TablellVl 



have been calibrated over a wide range of systems to yield 
a minimum of 4, 6 and 8 digits of relative accuracy in the 
total energy with LOOSE, GOOD and TIGHT respectively. 
Typically however, one additional digit is achieved for 
models that contain Hartree-Fock exchange. 



Accuracy 


Tmtrix 


T2E 


Tdist 


Thicu 


LOOSE 


10-" 


10" 




10-^ 


GOOD 


10-5 


10" 




10-5 


TIGHT 


10-® 


10" 


■10 


10-^ 



TABLE IV: Accuracy goals and corresponding thresholds that 
control precision of the MondoSCF linear scaling algorithms. 



